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Abstract 

The general tendency for species number (S) to increase with sampled area (A) 
constitutes one of the most robust empirical laws of ecology, quantified by species- 
area relationships (SAR). In many ecosystems, SAR curves display a power-law 
dependence, S oc A z . The exponent z is always less than one but shows significant 
variation in different ecosystems. We study the multitype voter model as one of the 
simplest models able to reproduce SAR similar to those observed in real ecosystems 
in terms of basic ecological processes such as birth, dispersal and speciation. Within 
the model, the species-area exponent z depends on the dimensionless speciation rate 
v, even though the detailed dependence is still matter of controversy. We present 
extensive numerical simulations in a broad range of speciation rates from v = 1CP 3 
down to v = 10~ n , where the model reproduces values of the exponent observed 
in nature. In particular, we show that the inverse of the species-area exponent 
linearly depends on the logarithm of v. Further, we compare the model outcomes 
with field data collected from previous studies, for which we separate the effect 
of the speciation rate from that of the different species lifespans. We find a good 
linear relationship between inverse exponents and logarithm of species lifespans. 
However, the slope sets bounds on the speciation rates that can hardly be justified 
on evolutionary basis, suggesting that additional effects should be taken into account 
to consistently interpret the observed exponents. 
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1 Introduction 



Species-area relationships (SAR) quantify ecosystem richness and, in partic- 
ular, the spatial variations of biodiversity. These curves measure the average 
number of species (S) present in a sample area (A) of a given ecosystem and 
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usual ly display a triphasic shape (iPrestonl . 119601 ; iRosenzweigl . 119951 ; IHubbelll . 
20011 ) . For small areas (below the dispersal range) and large areas (continen- 
tal scale), the number of species rapidly increases with the area; while for 
intermediate slower, sub-linear growth is observed. The intermedi- 

ate range is the most intriguing one and has gathered much attention since 



its discovery. Although many fun ctinoal forms have been 
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data in this intermediate regime (IHe and Legendre 
most common and widely accepted ones are the a lgebraic law S 
z < 1 and C a positive constant) proposed by lArrheniusI (jl92l|) 
Gleason ( 19221 )). and the l ogarithmic one S ~ Chi A due to Fisher et al 



(see also 



( 119431 ) . A recent survey by iDrakare et al.l (120061 ) , reconsidering most of the 



existing SAR studies from different ecosystems, shows that the former pro- 



pothesis cannot be a priori discarded, much efforts across the years (Preston. 


1962; 


MacArthur and Wilson, 


1967; Connor and McCov 




1979; 


Wriehtl. 1988: 


Kohn and Walsh. 1994: 


Durrett and Levin . 1996; Hubbell 


2001 


Chave et al.. 


2002; 


He and Legendre 


. 2002 


Martin and Goldenfeld, 


2006) have been de- 



voted to explain the observed values of the exponent z. Observations support 
the idea of a dependence of the exponent z on quantities s uch as latitude 



(lAUen and Gilloolyl . 120061 ) and body size of considered species ( IDrakare et al. 
20061 ). Notwithstanding observations and theoretical efforts, a satisfactory the- 
ory able to predict the value of the exponent in different ecological situations 
is still lacking. 



On the theoretical side, two distinct viewpoints on ecosystems organization 
correspond to different explanations for species-area relationships. According 
to the first, larger areas contain a lar ger variety of hab i tats a nd consequently 
can sustain a richer species diversity (IKohn and Walshl . Il994i ) . For the second 
viewpoint, species-area relationships are the outcome of demographic pro- 
cesses such as colonization, dispersal, speciation and local extinction, and do 
not need environmen tal diversity for their explanation (IMacArthur and Wilson! . 
19671 ; IHubbelll . l200ll ) . We should mention a third explanation, ascribing species- 
area relationships to sta tistical biases induced by the skewedness of species 
abundance distributions (IHe and Legendrd . 120021 ; iMartin and Goldenfeldl . l2006l ) 
As both niche-based and neutral-dispersal mechanisms are able to sustain di- 
versity, the hope is to extract information on t he importance of th e different 
classes of effects from the shape of SAR curves (jChave et al.l . |2002| ). 



We consider the voter model as the simplest prototype of neutral models able 
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to generate non-trivial species-area relationships (IDurrett and Levin I . Il996l : 
Hubbell l200ll ; IZillio et all 120051 : iRosindell and Cornell I2007T ). The model ac- 
counts, in a simple way, for the processes of birth, local dispersal and in- 
troduction of new species. Its main parameter, is, is a dimensionless number 
measuring the rate of appearance of new species — speciation events — in 
units of the death rate. The other ingredient is the dispersal kernel, which 
quantifies the probability for an individual of a species to colonize different 
locations in the ecosystem. Speciation and dispersal are enoug h to produce 



triphasic SAR curves resembling those observed in field data (IChave et al. 



20021 ; IRosindell and Cornell 120071 ). In particular, for local (short-range) dis- 
persal, the intermediate regime is well described by a power-law behavior 
S = CA Z , with an exponent z depending on the speciation rate v. We mention 
that the logarithmic function is reproduced by the voter model with global dis- 



persal, when indiv iduals can invade all loci of the ecosystem (j Coleman! 11981 



Chave et all |2002| ). The logarithmic law is also retrieved, for large dispersal, 
for areas being smaller than the dispersal range, indeed at these scales the 
dispersal appears as if it was long range. Understanding how the exponent z 
depends on the parameters of the model (in this case, mostly on is) is funda- 
mental to move a step toward the theoretical prediction of the variations of 
experimentally observed exponents in terms of ecological quantities. 



However, this dependence h as been source of s o me co ntroversy in the lit 



erature. In a seminal paper, IDurrett and Levin I (119961) proposed a formula 



accord ing to which, in the limit small is, z~l/ ln(l/ is). IRosindell and Cornell 
( 120070 sugge st a power l a w rela tionship between z and is. Finally, the scaling 
argument of IZillio et al.l (120051 ) predicts z approaching a finite value z ~ 0.2 
for vanishing is. These discrepancies have not yet been settled and, due to 
the weak dependence of z on v, a clean answer requires numerical simulations 
with v varying over several orders of magnitude. So far, only speciation rates 
v > 10 -6 were explored as simulations at lower (possibly more realistic) values 
of the speciation rate are computationally very expensive. 



In this paper we present results of simulations of the voter model with specia- 
tion rates varying in a wide range of values from v = 10~ 3 down to is — 10~ n , 
with the twofold aim of disentangling the low speciation rate behavior and 
examine an ecologically relevant range of parameters. Our findings are also 
useful to as sess whether ne utral predictions are consistent with realis tic speci- 
ation rates ( Hubbell l200ll ). a question which raised a heated debate (IHubbell 
20031 ; iRicklefsl . 120031 ). SAR curves resulting from our simulations are charac- 
terized by a power law behavior with exponent z, displaying a logarithmic 
dependence on the speciation rate and supporting de facto Durrett-Levin's 
scenario, even though with d ifferent numerical coefficients. In agreement with 
Rosindell and Cornelll (120071 ). we also found that the exponent z is essentially 
insensitive to the dispersal range implying that, accepting the hypothesis of 
the model, the observation of a species-area exponent imposes strong con- 
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straints on the rate of appearance of new species. In the Discussion section, 
we examine the plausibility of the model predictions on the basis of data avail- 
able in the literature. In particular, we consider the ^-values reported in the 
literature for different taxa and, due to the absence of reliable data on specia- 
tion times, we study how measured exponents depend on the lifespan, with the 
additional assumption that average speciation times and lifespans are linked 
by a scaling relation. The observed variations turn out to be much larger than 
those allowed in the framework of the model; we finally discuss which effects 
may be included to possibly achieve a quantitative description. 



2 Model 



We co nsider the voter model with mutation as defined by iDurrett and Levin 



( 119961 ). Individuals belonging to different species are placed at each site of a 
two-dimensional (L x L)-square lattice and evolve according to the following 
dynamics. At each time-step, a randomly chosen individual is killed, creating 
a gap which is immediately filled, with probability u, by an individual from a 
new species (not present in the ecosystem) — speciation event — or, with a 
probability 1 — u, by a new individual of an already existing (in the ecosystem) 
species chosen among those present in a neighborhood (that will be detailed 
below) of the site — birth/dispersal event. The dynamics is then advanced 
until the number of species in the ecosystem reaches a statistically steady 
value. Strictly speaking, the fact that empty locations ar e immediately colo- 



nized means that the birth rate is infinite (see discussion in Durrett and Levin 



( 119961 )). Therefore the basic time-step of the dynamics correspond to a death 
event, and thus the dimensionless parameter v represents the speciation rate 
a measured in unit of the death rate d. Equivalently, we can express v as 
the average species lifespan t divided by the average time between speciation 
events t^ s ' (we shall come back to this point in the Discussion section). 

As for the dispersal rule several options are possible. The simplest possibility 
is the nearest-neighbor rule, where the individual is replaced by one of the 
species present in the four neighbor sites with probability 1/4. We will refer 
to this in the following as the nearest-neighbor (NN) case. A more realistic 
choice is to use a generic dispersal kernel introducing the probability P(r) of 
a gap being filled by a species whose representative individual is at a distance 
r from it. We adopt a computationally simple instance by choosing the square 
kernel: we replace the individual with a copy of another individual randomly 
chosen in a square of side 2K + 1 centered on the gap. This choice does not 
represent a restriction as it has been shown that the relevant quantity is the 
averaged square dispersal distance and not the specific functional form of the 
kernel. For instance, a Gaussian and a square dispers al kernel with the same 



squared dispersal distance produce very similar SAR (IRosindell and Cornell! 
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2003). 



We stress that, independently of the dispersal rule, the model is completely 
neutral: all species (and individuals) undergo the same dynamics, as differences 
among species arise only due to demographic stochasticity. 

Simulations have b een efficiently perfo r med by using the dual representation 
of the voter model (iHoiley and Liggettl . Il975l ). providing a way to reconstruct 
the asymptotic configuration of the ecosystem by tracing backward in time its 
evolution. An important advantage of the dual representation is that it recon- 
structs the genealogy of each individual up to the speciation event originating 
its species, meaning that the system is ensured to have reached equilibrium. 
Moreover, it allows to implement open boundary conditions: the genealogy of 
an individual can be reconstructed also when its ancestors are outside the sim- 
ulated are a, which can thus be c o nside re d as a sample o f a vir tually infinite 
ecosystem iRosindell and Cornell! (120071 ): iRosindell et al.l (120081 ). This means 
that we can interpret v as a bona fide speciation rate, since immigration from 
outside the system is included in the birth-dispersal process. However, long- 
range immigration events qualitatively different from local dispersal (i.e. seeds 
transported by birds) can be modeled as an higher "effective" speciation rate 
v. We managed to optimize the algorithm to simulate the model for very low 
speciation rates, down to v = 10~ n . Details on how the simulations have been 
perform ed and the stati s tics h ave been collected can be found in Appendix |A| 
see also IRosindell et al.l (120081 ) for other possible improvements of the coales- 
cence algorithm. As for the dispersal, we explored both the NN and the square 
kernels, for the latter K has been varied in the range K = 3 — 64, though we 
shall mostly present the results for K = 7 (see the discussion in the next 
section). 



3 Numerical results 



We begin studying SAR curves obtained at fixed dispersal range (K = 7) 
and varying v, as shown in Fig. [TJ All curves display a fast growth for small 
areas with a crossover, for areas of the order of the dispersal kernel (A 
K 2 ),to the power-law regime. The final regime where the number of species 
becomes linear with the area can be detected only for rather large speciation- 
rate values, 10~ 5 < v < 10 -3 ; to observe it at lower values of v much larger 
simulation samples would be required. In the inset, we plot the "local species- 
area exponent" for each curve, d(ln S) / d(ln A) , which clearly shows that the 
smaller v the smaller the exponent becomes and the larger is the range of 
scales where a well defined power-law behavior establishes. Finally, when the 
parameter v is not too small, it is possible to o bserve also the final linear r egime 



which occurs for areas much larger than v 1 (IDurrett and Levin I . Il996l ). 
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Fig. 1. Dependence of SAR curves on the speciation rate: S vs A for different spe- 
ciation rates v G [10 -9 : 10~ 3 ] obtained with the square kernel with K = 7 in 
a simulation sample of side L = 2048, and averaging over 100 independent real- 
izations. Note the triphasic shape clearly observable for larger values of v. Inset: 
Logarithmic derivatives of the SAR curves, d(lniS)/d(ln A) vs A, a plateau identifies 
the intermediate regime and the plateau value the exponent z. Note that z increases 
with v and the intermediate regime enlarges in width at decreasing v up to invading 
almost all the simulation sample for small v values. 

Figure [2] (left and middle panel) exemplifies the behavior of species-area curves 
at fixed v and different dispersal range K. At increasing the dispersal range 
the onset of the power law regime shifts at larger areas, apparently with- 
out affecting the exponent. A more careful analysis of the local exponents 
d(lnS')/d(ln A), shown in the right panel, detects a dependence of the value 
of the exponent on the dispersal range when this is small, K < 5, including 
the NN case. 



On the other hand, when K > 5, we did not observe any appreciable correc- 
tions to the value of the exponent. The indepen dence against variations of K , 
when it is large enough, has been quantified by lRosindell and Cornell! (120071 ). 
who have shown that curves obtained with different (not too small) K can be 
rescaled on a universal function of A and v only via the transformation: 



S = f(A,u,L)=K r <f ) (A/K r 



v) 



characterized by the scaling exponent r w 1.97. We checked that this rela- 
tion holds also with the small values of v that we studied, for instance the 
insets of Figure [2] (left and middle panel) show it for v = 10~ 5 and v = 10 -8 . 
We will then study in the followin g the NN and the K = 7 cases, the for- 
mer being that originally studied by iDurrett and Levin I (119961 ) and the latter 
being representative of the behavior of the model for large average dispersal 
distances. 



We now turn to the main results of this paper about the dependence of z 
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Fig. 2. Dependence of SAR curves on the dispersal length. Left and Middle panels: 
S vs A, for v = 10 -5 and v = 10 -8 obtained with the square kernel with several 
values of K in a simulation sample with L = 2048. Note that in the intermediate 
regime, which starts and ends at different area values by changing K, the slope of 
the curves is the same indicating that z is independent of K. Insets show the curves 
rescaled with the transformation ([I]). Right panel: local slopes d(ln S) / d(ln A) vs A 
for v = 10 6 and different kernels: NN and square with K = 1,3,5,7,11. Note the 
tiny dependence of the exponent (the region between the two horizontal lines) on 
the dispersal length. For K > 5 the plateau region (lower horizontal line) does not 
change anymore apart from an horizontal shifting of its onset towards larger areas. 

on v. In Figure [3] (left) we show the exponent function of the speci- 

ation rate v (see Appendix [S] for a discussion on how we estimated z). We 



observe a clear discrepancy with previous predictions ( IDurrett and Levin 



1996 : Rosindell and Cornell! 2007 ) (also shown in the picture). In particular, 
for v <C 1, we found the data to fall into a straight line when plotting 1/z vs. 
ln(z/) (Fig. [3] right), suggesting the following functional dependence 



q + m ln(z/) 



(2) 



by which we obtained a best fit to the data with q « —3.3 and m ~ —0.72. 
To compare our results with p revious studies of these model s, notice that the 
power-law fit suggested in by iRosindell and Cornell! (120071 ) agrees with the 
data in the same range of speciation rate values, i.e. u > 10 -5 . Deviation 
from a power law behavior are clearly observed for lower values of v, where 
the data also rule out the saturation at z « 0.2 pr e dicted by IZillio et al. 



( 20051 ). Actually, our fit confirms IDurrett and Levin ] (Il996h prediction of a 
logarithmic decay of z with u, up to corrections order (9(ln(ln(z/))). However, 
the fitting parameters m and q for both th e square kernel with K = 7 and the 



NN kernel are very different from those of IDurrett and Levin I (see caption of 



Fig. [3]). We conjecture that th e differences in prefac tors could be caused by two 
different assumptions used by IDurrett and Levin I to derive the dependence of 
z on v. The first is about pre-asymptotic effects: the statistical results used by 
Durrett and Levin are strictly valid only when t — > oo which requires v — > 0, 
while finite-time corrections may affect the exponent value. In this respect, 
also for our data the v — > limit seems to be crucial for the validity of the 
fit 0. The second is the assumption that a power law regime establishes from 
A = 1 to A = v~ x . Conversely, we observe the onset of the power law for areas 
being slightly larger than 1 even in the NN case. Moreover, the crossover to 



7 




Fig. 3. Dependence of the species-area exponent on the speciation rate, (left) SAR 
exponent z versus the speciation rate v for the NN and square kernel wit h K = 7. 



For comparison it is also shown the prediction by Durrett and Levin 1 ( 19961 ) and the 
power-law dependence suggested bv lRosindell and Cornell (|2007l ). (right) Same data 



as left but shown plotting 1/z vs ln(z/). The straight lines reports the best fit with 
© with parameters: m = -0.48 ± 0.02, q = -1.4 for NN data, m = -0.72 ± 0.02, 
q = —3.2 for square kernel data (K = 7) and m = —0.127 ± 0.002, q = 2.10 for 
Durrett-Levin tabulated data. Errors on the exponents are of the order of symbol 
sizes, see Appendix A for details. 

the linear asymptotic regime begins for areas quite smaller than 1/z/. 

It should be noticed that discrepancies in the numerical factors have profound 
implications when the model is used to estimate a speciation rate from an 
observed species area exponent. The logarithmic dependence of z on v makes, 
in fact, v exponentially dependent on z. We will discuss in the next section 
how this dependence ca n be compared with experimental data. It should also 



be remarked that both iDurrett and Levin I prediction and equation (j2J) are 



valid for small values of v and can lead to incorrect results, such as negative 
z, for v close to 1. 

All simulations so far presented have been performed with open boundary con- 
ditions, which are appropriate when the sample is a homogeneous portion of a 
much larger ecosystem. However, closed boundary conditions can be of interest 
for modeling confined ecosystems such as islands. Intuitively, open boundaries 
allow new species to immigrate into the sampled system from the external 
infinite ecosystem, independently of the speciation events. Closed boundaries 
exclude this possibility and are thus expected to reduce the exponent z and, 
in general, species richness. Fixing the speciation rate v the decreasing of z 
becomes more and more efficient as the system size decreases, and at fixed 
size the effect is the stronger the smaller is v. The closed boundaries effects 
becomes dramatic for simulation samples with A ~ 1/z/, and, in particular, we 
observed that for v such that A < 1/z/, biodiversity is definitely lost, i.e. z = 
(see FigH]). Notice that when the system size is large (i.e. A > 2048 2 ) and the 
speciation rate is not too small (i.e. v < 10~ 4 ) the exponent is essentially the 
same of the open boundary case. We also remind that, if islands are modeled 
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Fig. 4. z vs v for the square kernel with K = 7 with closed boundary conditions 
(CBC) at L = 512 and L = 2048 compared with the result with open boundary 
conditions (OBC) and L = 2048. Note how the decrease of z with v becomes more 
and more dramatic at decreasing the size. 

by closed boundaries, the par ameter v should be meant to include the im- 
migration rate of new species (IMacArthur and Wilson! 119671 ). since dispersal 
from outside the system is forbidden in this case. 



We stress here that the implementation of closed boundary is a very simplis- 
tic way of describing confined ecosystems and more sophisticated effects may 
be relevant in these cases. For example, it has been shown that the dynam- 
i cs at the edge of eco systems can be quite different from that in the bulk 



( iLaurance et all 120061 ). This means that "open" situations, when the sam- 



ple is part of a much larger and homogeneous system, provide a much safer 
comparison between models and field data. 



4 Comparison with empirical data and discussion 



Species- Area relationships have been subject of intense experimental research 
in a variety of ecosystems, and the range of variability of t he exponent describ- 



ing th e intermediate regime goes from z ~ 0.05 in bac teria (IHorner-Devine et al. 



20041 ) to z = 0.4 — 0.5 in some plants community (see lDrakare et al.l (120061 ) for 
an exhaustive review of field observations). According to our results, the voter 
model with speciation displays a variability of z in the same range when the 
speciation rate is allowed to vary over several orders of magnitude. It is thus 
tempting to go in the direction of a more quantitative comparison between 
field data and the model results. 

As a first step, we assume that the model is able to describe the main fea- 
tures of groups of trophically similar species and explore the consequences of 
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this assumption. This requires that a speciation rate and a dispersal range 
for the whole group o f species can be properly defin ed, although we know 
that speciation rates (jMariakeva and Gorshkovl. I2004T) and d ispersal ranges 
flNathan and Muller-Landaul . 120001 " iKinlan and Gained . 120031 ) may have sig- 
nificant variations from species to species. However, in the model the exponent 
is essentially independent on the dispersal range and depends only logarithmi- 
cally from the speciation rate, so that these variations might be disregarded 
treating all species in the group as having the same "average" dispersal range 
and speciation rate. 



As far as dispersal is concerned, we only found a tiny dependence for very short 
dispersal range, around K < 5. Above these values th e exponent is indepen- 
dent o f the dispersal range confirming pre vious results iRosindell and Cornell 
( 200 73) • Realistic average dispersal ranges (INathan and Muller-Landaul . |2000| ; 
Kinlan and Gained . 120031 ) are certainly far from the short dispersal case, due 
to animal motility or wind for seeds. Therefore, we assume that the dispersal 
range of real groups of species is always in the range where the exponent is 
dispersal-independent. It is however worth remarking that the dispersal range 
can still affect the spatial biodiversity via the power-law prefactor, whose in- 
crease can lead to a large number of species that, when z is small, increases 
very slowly with the area. In this respect, the model outcomes are in contrast 
with interpretations of low v alues for z in bacteri a as an effect of large dis- 



persal distances as argued in iDrakare et al] (120061 ) and iHorner-Devine et al. 
(120041 ) . 



What about speciation? Unfortunately, we do not have ecological data al- 
lowing us to directly estimate the frequency of speciation events. Data from 
fossils su ggest an a verag e speciation rate on Earth of about three specie 
per year (jSepkoskil . Il998l ). but it is hard to infer from this number 
sonable rate for a living system . Also estimates based on mutation rates 
( jMariakeva and Gorshkovl . 120041 ) could be flaw ed due to genetic bottlenecks 
and phenomena like horizontal gene transfer ( Jain et al. , 1999). Moreover, 
as discussed in the Model section, the parameter v should be interpreted as 
an "effective" speciation rate, incorporating also long-range dispersal events. 
Within the model framework, our results show that species-area exponent 
and dimensionless speciation rate v are related even when the latter is very 
small, implying that an observed value of the exponent z would predict the 
rate of introduction of new species v. Remarkably, the existence of positive 
correlations of these two quantities is consistent with observational results. 
As an example, it is k nown that close to the equator species-area expo- 
nents tend to increase (iDrakare et all 120061 ) toget her with speci ation rates 
( Allen and Gillooly . 20061 ) and overall biodiversity ( Stevens . 1986 ). 



In order to test the ecological plausibility of the relation between z and u, 
we make use of the definition of v as the ratio between the speciation rate a 
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and the death rate d. From Eq. ([2]) and separating the contribution from the 
variation in the speciation rate from that of the variation in the death rate, 
we have 

- = q + mln(V) = q + m[ln(<r) — ln(cf)] , (3) 

where the arguments of the logarithm are made dimensionless by measuring 
them in the same units. To ease the interpretation, we recast this equality in 
the time domain using the lifespans t = 1/d and the average time between 
speciation events f( s > = 1/er: 

i = g + mpn(*)-h(* W )]> ( 4 ) 

The first term on the right hand side accounts for the variation in z due 
to the lifespan which is, of course, much easier to estimate than the term 
due to speciation time and can still be important and informative. Indeed 
there are evidences that taxa having a shorter gene r ation time have generally 



lower species area exponents (jHorner-Devine et all 12004 ; iGreen et all 12004 



Zhou et al. . 20081 ) (we recall that m is negative). We thus study how the 



inverse exponent 1/z varies with the logarithm of the lifespan. The results of 
this analysis are presented in figure for data obtained from the literature 
(see Appendix [B] for a description of how data have been collected), which 
shows that a linear relationship fits rather well the data, with an observed 
slope m meas = —1.76 ± 0.13 (dashed line in the figure) which is different from 
m pa —0.72 predicted by the voter model. 

The fact that for species-area exponents measured in field data we found 
1/z oc ln(t) suggests a scaling relationship between speciation time and lifes- 
pan, i.e. that ~ P, so that m meas = m(l — 7), as clear by substitution in 
the previous formula. We do not have any a priori explanation for justifying 
a power law dependence of the speciation time on the lifespan, apart from 
the observation that the variations of many ecologi c ally r elevant rates among 



species are governed by scaling laws (IBrown et all |2004|). We are not aware 



of independent estimation of the dependence of the speciation time on the 
lifespan so to confirm or reject the outcome of our analysis. 

However, the relation m meas = m(l — 7), with m and m meas fixed by the 
voter model and field data respectively, yields a negative 7. This result is in 
contrast with biological expectations as it would imply, e.g., a speciation time 
for bacteria much longer than the one for trees, which is hard to justify bio- 
logically. Reasonable expectations would have been < 7 < 1. The limiting 
case of 7 = 1 is the trivial case in which speciation time is proportional to the 
lifespan. This would have lead to m meas = 0, i.e. same v = t/t^ and z for all 
taxa. The other limiting case is 7 = 0, which is plausible when the possibility 
of creating a new species is triggered by some external mechanism, like the 
availability of new niches, which is not strongly correlated to any particular 
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Fig. 5. 1/z vs logarithm of lifespans (measured in days) for several taxa. Dashed 
line is the best fit, giving a slope m meas = -1.76 ± 0.13 (reduced x 2 = 1-04). 

feature of the species. Another justification could come from co-evolutionary 
mechanisms: species having very different lifes pans can still evolve on simila r 
timescales due to their ecological interactions ([Thompson and Jeremyl . Il992h . 
Actua lly, co-speciation is k nown to occur in some cases of host-parasite sys- 
tems ( Clayton et all 2003 ). In this case, one would find m meas = m. 



The inconsistent value of 7 we obtained can be interpreted either as a failure 
of the basic assumptions of the neutral model and thus of its inadequacy in 
describing realistic ecosystems or as the need of additional mechanisms to be 
included in the framework of dispersal models. The quality of the linear fit 
shown in Fig. [5] would suggest to opt for the second interpretation, even if 
the linear dependence of 1/z on log(t) needs to be tested with further mea- 
sures. In recent years, several attempts of relaxing the strong assumptions of 
neutral models have been tried. The results of these models are pretty ro- 
bust with respect to modifications of so me hypothesis such as the saturation 
of the resources (jEtienne et all l2007bl ). In spatial models, it has been also 
observed that the introduction of trade-offs does not have a dramatic effect 



on species-area exponents (IChave et al.l . 120021 ) . Therefore, it is reasonable to 
search for other elements in the model which can lead to a failure in repro- 
ducing the observed data. In particular, the assumptions of a point speciation 
mode (i.e. the fact that each individual has a fixed probability to speciate) is 
known to be crucial and the results may change dramatic ally when considering 
"fission modes " , corresponding to al lopatric speciation (iHubbelll . l200ll . 120031 ; 



Ricklefd . 120031 ; lEtienne et al.l . l2007al ). Indeed, it seems like the most impor- 
tant assumption of neutral model s is that new species ent er the system with a 
population of a single individual (jZillio and Conditl . 120071 ) . This could explain 
why speciation rates predicted by neutral models with point speciation may 
look unrealistically high: new species are introduced with only one individual 



12 



and they have an high probability of going extinct before being able to grow. 
In other words, there could be a discrepancy between the parameter v in the 
model and the experimentally observed speciation rate. The effect of realistic 
speciation mechanisms on neutral and more general dispersal models could 
be key to understand the puzzle of the observed variation of the exponents 
among different taxa. 

In conclusion, simulations of the multi-type voter model for low values of 
the speciation rate show a clear logarithmic functional dependence of the 
specie-area exponent on the speciation rate and independence on the dispersal 
kernel (provided it is not too short ranged). Analysis of field data support a 
logarithmic dependence of the exponent z on the timescales of the problem, 
though with a prefactor which is incompatible with that found in the model. 
Our analysis points out that more refined models should allow larger variations 
in the exponent z in order to be consistent with observational data. 
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A Simulations and data analysis 

By means of the dual representation of the voter model, the model becomes 
equivalent to a system of coalescing random walkers, where a birth/dispersal 
event corresponds to collision-coalescence of two walkers and a speciation event 
to the death of a walker. Therefore, as time proceeds, less and less surviving 
walkers should be accounted, speeding up the simulation. This allowed us to 
simulate lattices L x L with L = 1024, 2048 and L = 4096 for v e [10~ n : 
10 -3 ]. In order to test the effect of the system size when boundaries are present, 
we have also used L = 512; in this case we generalized the algorithm by just 
refusing all moves causing the exit of a walker from the simulation domain, 
thus constraining the walkers to remain inside the initial grid. To embank the 
simulation bottleneck due to the initial presence of L 2 walkers we optimized 
the walker collision detection by means of a look-up table. For each v and L 
we repeated the simulation many times with different seeds for the random 
number generator, typically from 100 to 150 - 300 for L = 2048, 4096 and 
L = 1024, 512, respectively. For the lower values of v simulations get very 
slow and a lower number of realization was used, typically from 20 to 60 for 
v < 10~ 10 . Once the species occupancy patterns are obtained SAR curves 
are derived by averaging the number of species in non-overlapping squares of 
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side y/A = 1, . . . , L whose union completely covers the simulation grid. So 
that averages are performed both over the number of sampled areas in each 
realization and over different realizations. Statistical errors on the average are 
also computed. The exponent z characterizing the power law growth of S 
with A is then estimated fitting by a linear least square method the function 
In S = q + z In A with A e [A min : A max ] chosen at the beginning and the end 
of the intermediate regime, respectively. The fit were performed by minimizing 
the reduced x 2 function (i.e. normalizing the \ 2 with the number of degrees 
of freedom) but constraining the minimal number of points to be considered 
(from 5 to 15 depending on the extension of the intermediate range). As the 
least square error is smaller than the variability of the fitted z at changing 
the minimal number of constrained points, we set the error on the estimate 
as such variability. In Fig. [3] errors are comparable with the symbol size. 

The quality of the fit is then compared by a direct inspection of the local 
slopes (logarithmic derivatives, i.e. dhiS/dlnA vs In A) of the SAR curves. 



B Details on observational data 



Da ta presented in Fig. [ 5 ] are b ased on the collection of exponents presented 



by iHorner-Devine et al.l (120041 ) . We avoided presenting the z- value for plants 
since it varies a lot among d ifferent studies (see supplementary information of 
Horner- De vine et al. ( 20041 )). As a representative of long-lived organisms we 



have chosen tr opical forest tre es which are well studied and w e assumed for 
them z = 0.27 riLonsdakd . [l~999h and average lifespan « 60 years (ICondit et al.l 



19991 ). The exponent z values for butterflies, earthworms, birds and ciliates are 
the same of the original refere nce; the value for bact e ria z = 0.05 is an average 
between the value z = 0.04 in Ho rner-Deyine et al.l (120041 ) and the value z = 
0.06 in a more recent study by Zhou et al. ( 20081 ). As far as the other lifespans 
are concerned, we must stress that it can vary much from species to species 
and, in most cases, it is hard to find in the literature good statistical studies. 
However, due to the logarithmic dependence, the fit is not be so sensitive 
to errors in the estimates as far as the order of magnit ude is correct. W e 
assumed an average lifespan of 2 years for birds (see, e.g., ISpeakmanl (120051 )). 



Despite their short average lifesp an in their adult stage, b utterflies usually 
have a few generations per year (IGilbert and Singerl . Il975l ); we assumed an 
average of two, correspondi ng to t = 0.5 years . Othe r estimated lifespans 
are 1 year for earthworms (Edwards and Bohlenl . 119770. 10 days for ciliates 
(IJensen and Verschoorl . I2004J ) and 2 days for bacteria (jClarholm and Rosswalll . 



19801 ). 
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